Power law relaxation and glassy dynamics in Lebwohl-Lasher model near 

isotropic-nematic phase transition 
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Orientational dynamics in a liquid crystalline system near the isotropic-nematic (I-N) phase tran- 
sition is studied using Molecular Dynamics simulations of the well-known Lebwohl-Lasher (LL) 
model. As the I-N transition temperature is approached from the isotropic side, we find that the 
decay of the orientational time correlation functions (OTCF) slows down noticeably, giving rise to 
a power law decay at intermediate timescales. The angular velocity time correlation function also 
exhibits a rather pronounced power law decay near the I-N boundary. In the mean squared angular 
displacement at comparable timescales, we observe the emergence of a subdiffusive regime which is 
followed by a superdiffusive regime before the onset of the long-time diffusive behavior. We observe 
signature of dynamical heterogeneity through pronounced non-Gaussian behavior in orientational 
motion particularly at lower temperatures. This behavior closely resembles what is usually observed 
in supercooled liquids. We obtain the free energy as a function of orientational order parameter by 
the use of transition matrix Monte Carlo method. The free energy surface is flat for the system 
considered here and the barrier between isotropic and nematic phases is vanishingly small for this 
weakly first-order phase transition, hence allowing large scale, collective and correlated orientational 
density fluctuations. This might be responsible for the observed power law decay of the OTCFs. 
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I. INTRODUCTION 



Liquid crystalline systems often exhibit interesting dy- 
namics apart from the rich phase behavior. Surprisingly, 
dynamics of such systems have traditionally been probed 
only at rather long timescales (nanoseconds to millisec- 
onds) P, 0, • Recently, Fayer and coworkers have in- 
vestigated the dynamics in the isotropic phase of ther- 
motropic liquid crystals over a wide range of timescales 
using optical Kerr effect (OKE) measurements 0,|EEI0- 
At short to intermediate timescales, they have observed 
a pronounced power law decay of the time-dependent 
OKE signal near the I-N phase transition. At the inter- 
mediate times (several nanoseconds) the decay becomes 
even slower, almost appearing as a plateau on a log-log 
scale Q . The exponential decay predicted by Landau-de 
Gennes theory is observed only at the longest timescale 
(> 10 ns). 

Subsequent molecular dynamics simulations of a 
calamitic system (comprising of rod-like molecules) with 
the Gay-Berne pair potential were found to reproduce the 
power law decay of orientational time correlation func- 
tions (OTCF) PI. From a very recent computational 
study, which deals with a calamitic system, a discotic 
system (comprising of disc-like molecules) and to a lim- 
ited extent the Lebwohl-Lasher lattice model, it appears 
that this power law decay at short-to-intermediate times 
might be a rather general phenomenon in thermotropic 
liquid crystals 0- This observation has gained support 
from a recent finding of power law decay at relatively 



short times in an idealized calamitic liquid crystal model 
with length-to- width ratio 5-6 ^3] • It has been observed 
that many aspects of the orientational relaxation behav- 
ior outlined above bear close resemblance to what is ob- 
served in supercooled molecular liquids near the glass 
transition temperature fH 1^. ITU]. In particular, the de- 
scription of orientational relaxation with a power law de- 
cay at short times and exponential decay at long times 
is strikingly similar 7]. However, the origin of such a 
rich dynamical behavior may be quite different in the two 
cases • It would be especially interesting to explore the 
free energy surface with respect to orientational density 
fluctuations in search of the origin of the slow dynamics 
near the I-N transition. 

Comprehensive understanding of this rather exotic dy- 
namics of thermotropic liquid crystals spread over al- 
most five decades of time is a challenging task, particu- 
larly because of the anisotropic nature of the interaction, 
which makes the theoretical analysis difficult. Computa- 
tional approaches have become very much useful in this 
regard as we gain control over the microscopic interac- 
tions and try to understand their manifestation into the 
macroscopic behavior. In this spirit, this communication 
attempts to continue the investigation of the Lebwohl- 
Lasher model. 

Lebwohl-Lasher (LL) model [lj is essentially the lat- 
tice version of the Maier-Saupe model. In this model, the 
molecules being fixed on a simple cubic lattice lose their 
translational degrees of freedom and can only rotate. The 
total interaction energy is given as follows: 
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strength of interaction between the rotors. It has a fixed 
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value e for the nearest neighbors and otherwise. 

LL model has been rather popular since its inception in 
computer simulation studies of liquid crystals. Its merit 
lies in the inherent simplicity with which it clearly estab- 
lishes some of the very essential features of orientational 
ordering in liquid crystalline systems 0, 0, 0, 0, 0] . 
Simple form of the anisotropic interaction and the ab- 
sence of translational degrees of freedom makes the sys- 
tem particularly easy to study. Most of the earlier stud- 
ies involving LL model used Monte Carlo (MC) methods, 
which can not predict the real dynamics of the system. 
In view of this, here we have undertaken molecular dy- 
namics (MD) simulations to investigate the dynamics in 
the LL model. 

Our MD simulations of the LL model show that as the 
transition temperature is approached from the isotropic 
side, the decay of the orientational time correlation func- 
tions slows down noticeably, giving rise to power law de- 
cay at intermediate timescales. Another interesting re- 
sult is the emergence of subdiffusive orientational motion 
essentially in the same temporal window where the power 
law decay sets in. In addition, the sub-diffusive motion 
is followed by a superdiffusive regime before the long- 
time diffusive behavior sets in. Moreover, we have used 
the recently developed transition matrix Monte Carlo 
(TMMC) method to obtain the free energy profile for 
the system as a function of orientational order parame- 
ter. We find the barrier to be vanishingly small for this 
weakly first-order transition allowing large fluctuations 
in the orientational order. 

The organization of the rest of the paper is as follows. 
In the next section, we discuss the details of simulations. 
In section Hill we present the results on orientational re- 
laxation. Section Hvl deals with the TMMC method and 
the results it yields on the free energy calculation. Sec- 
tion^concludes with summary and a few pertinent com- 
ments. 



II. SIMULATION DETAILS 

The system under study consists of a (10 x 10 x 10) 
simple cubic lattice with one rotor fixed on every lat- 
tice point. Only the nearest neighbors interact through 
the orientation dependent potential given by Eq. 
Throughout we have used dimensionless temperature T* 
and time t* (for the MD studies) defined by |l3j : 

T* = k B T/e 

t* = tie/I^) 1 ' 2 (2) 

Here both the parameters e and I± (moment of inertia 
with respect to the axis perpendicular to the molecular 
axis) have been taken to be unity. 

We have performed Molecular Dynamics (MD) simu- 
lations in microcanonical (NVE) ensemble using velocity 
Verlet algorithm. We have used time step of St* = 0.002 



in reduced unit and have obtained energy conservation 
upto fifth place of decimal throughout the simulation. 
We have scaled the velocities at every 100 steps for ini- 
tial 10 5 steps to equilibrate the system at a particular 
desired temperature and stored the trajectory for analy- 
sis after allowing the system to evolve without scaling for 
another 10 5 steps. The standard deviation in tempera- 
ture during data acquisition was of the order of 0.01 for 
all temperatures. Periodic Boundary Conditions (PBC) 
have been applied to remove the surface effects. Earlier 
MC studies showed pronounced system size dependence 
for this model. The same applies for the MD simulation 
also. The size effect is particularly important near the 
transition temperature as the correlation length tend to 
diverge. The transition becomes sharper with increas- 
ing system size. The transition temperature {Tin) and 
the free energy barrier follow certain finite size scaling 
laws 0| . We have found the dynamical features reported 
here to be qualitatively similar for different system sizes. 

III. RESULTS 

The phase behavior for LL model has been well ex- 
plored using Monte Carlo (MC) methods for pretty large 
systems. But since there are relatively less number of 
Molecular Dynamics (MD) studies, we have computed 
the average order parameter at various temperatures us- 
ing MD trajectories. The value of orientational order pa- 
rameter is determined by diagonalization of the ordering 
matrix O 



Qa(3 = Wjy ^ [3^0,6,0 - 5 a p], (3) 
i=l 

where ei a is the a-th Cartesian coordinate of the unit 
vector (ej) specifying the orientation of the i-th molecule. 
The largest eigenvalue and corresponding eigenvector 
give the orientational order parameter and the director 
respectively for a particular configuration. The thermo- 
dynamic order parameter {S) is obtained by averaging 
over the simulation trajectory. As shown in Fig. ^ the 
data obtained from the MD simulation overlaps well with 
the result obtained from our MC simulation. 

Note that MD simulation gives us the value of the tran- 
sition temperature to be ~ 1.14. This differs from the 
more accurate value (~ 1.1232) obtained from MC sim- 
ulations 0, a t the second place of decimal. This 
is acceptable because of the fluctuation in temperature 
present in microcanonical (NVE) MD simulation and the 
standard deviation of the order of 0.01 at all temper- 
atures. The fluctuation in temperature is particularly 
high near Tin- We have used the value obtained from 
our MD simulation to explain our other observations. 

Often the second rank order parameter is not enough 
to decide if true long range order exists in the system. For 
that purpose, orientational correlation functions (G;(r)) 
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FIG. 1: The average orientational order parameter indicating 
the I-N transition at T* = 1.14. Empty circles indicate the 
MC results and filled circles indicate the MD results. 
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FIG. 2: The nature of spatial decay of the two particle orien- 
tational correlation function G2(r) clearly demonstrates the 
existence of long range order in nematic phase. The correla- 
tion function decays to very small value as isotropic phase is 
reached. 



are particularly useful. The set of correlations can be 
defined as expansion coefficients of the rotationally in- 
variant pair distribution [T?) : 

or i 1 

G(r,/3 12 ) = G°°(r)^-^G i (r)P i (cos%), (4) 

where Gg°(r) is the particle center distribution, (9y is 
the angle between rotors % and j, and Pi is the Legendre 
polynomial of rank I. For a simple cubic lattice, 



where p is the density and Z], is the number of neigh- 
bors at rfc. Hence, Gi(r) = (Pi(cos8ij)) r gives orienta- 
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FIG. 3: Demonstration of power law regimes in the log-log 
plot of single particle orientational correlation function C| (t) : 
Fitted power law regimes are shown in thick lines. The decay 
is almost exponential for the highest two temperatures, i.e. 
the power law component becomes very small. 



tional correlation between two rotors separated by dis- 
tance r. Figure [21 demonstrates that if true long range 
order is present, G2(r) decays to a plateau with value 

(p 2 ) m 



A. Power law decay in OTCFs 

Single particle OTCF gives a temporal measure of the 
loss of the memory of a single particle of its own ori- 
entation in the environment created by the surrounding 
molecules. The single particle OTCF of rank I is defined 
as: 



cm = 



(Eifl(ei(O)-ei(t))) 
<Eiflte(0)-e<(0))) 



(6) 



where e,-(t) is the unit vector denoting the orientation of 
i-th molecule at time t. Since LL model has up-down 
symmetry, Cf (t) would be physically meaningful. 

In Fig. |3| the log- log plot of Cf (£) versus time is shown 
at different temperatures approaching Tin- An inter- 
esting step-like decay is clearly evident. The power law 
behavior emerges just above the transition temperature, 
i.e. in the isotropic phase close to the I-N transition and 
continues in the nematic phase. It has been rather dif- 
ficult to fit the entire decay by a single function due to 
the oscillatory regime in the nematic phase and rather 
dominant plateau appearing just before the power law 
regime. Hence, the OTCF has been fitted only beyond 
the initial plateau using the following functional form. 



C 2 s (i) = A + cxp(-t/T)(t/r p )- a 



(7) 



Here, the parameter r gives the timescale correspond- 
ing to the long time exponential decay, the parameter t p 
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signifies the timescale of the power law decay and the 
parameter a gives the power law exponent. The values 
for r, t p and a are shown in the following table. Note 
that the nature of the decay of the single-particle OTCF 
is strikingly similar to the experimental results obtained 
by OKE measurements. 
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Tp 
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1.13 


66.17 


0.002 


0.33 


1.15 


12.86 


0.1 


0.82 


1.16 


3.67 


0.07 


0.64 


1.18 


1.69 


~ 0.0 


0.001 


1.21 


0.84 


~ 0.0 


~ 0.0 



We have an interesting observation that the timescale 
of the exponential decay at the longer times becomes pro- 
gressively larger with decreasing temperature, whereas 
the values of t p clearly indicate that power law decay 
exists only near I-N phase boundary and this becomes 
transient as one moves away from the transition region. 
The values of the exponent a follows a similar trend. 
The OTCFs for the highest two temperatures are almost 
purely exponential beyond the plateau. We interpret the 
above result in the following fashion. As the temperature 
approaches Tin from the isotropic side, the single par- 
ticle orientational relaxation becomes slower due to the 
strong coupling of rotational motion with the surround- 
ing molecules, which result in an orientational caging ef- 
fect of the rotors. The formation of orientational cage 
or pseudo-nematic domains (within the isotropic phase) 
surrounding a molecule, gives rise to the retention of the 
memory of the previous orientation for a longer time, 
resulting in the slow relaxation process. Thus, the col- 
lective relaxation becomes important near the transition 
point. 

Mode coupling theory (MCT) can be used to ob- 
tain a semi-quantitative theoretical description and also 
physical insight of the slow dynamics near phase tran- 
sitions O 120, EH- MCT starts with splitting the fre- 
quency (z) dependent memory function, here the rota- 
tional friction, into two parts: 



T R (z) = T RB (z) + T Rp (z) 



(8) 



where T R b(z) is the short range or binary part which 
decays on a short timescale and usually originates from 
collisions between the molecules. Whereas, T Rp (z) is the 
collective part, deriving contribution from collective cor- 
relations. 

In the present lattice model there is no collisional con- 
tribution to the rotational friction. Therefore, T R b(z) 
can be set to zero and this is in fact the reason for the 
large contribution of the initial inertial decay. Friction 
is small at short times. However, T Rp (z) exhibits singu- 
lar features at small frequencies due to the emergence of 
long range orientational correlations near the I-N phase 
boundary. It was shown elsewhere that at low frequency, 



-0.8 



C3 

d 

.§> .12 





* 

- T = 1.15 


1- 


T = 1.18 



-1.4 



01) 

o 



-1.6 



-1.8 



0.4 



0.6 



log 



10 



0.8 
(t) 



1.2 



FIG. 4: Power law is exhibited in the derivative of collective 
OTCF, which is essentially the OKE signal. Fitted power law 
regimes are shown by thick lines. 



the frequency dependent friction develops a rapid growth 
which can be represented as Q: 



r<» ~ a z - 



(9) 



Mean field treatment gives a = 0.5. In general, invoking 
the rank (I) dependence of the memory function, the sin- 
gle particle OTCF can be expressed as [2I IH ElEHj : 



z + 



l{l + l)k B T 

J(* + r,(*))J 



(10) 



As discussed by Gottke et al, this expression can be 
Laplace inverted to obtain a power law decay in C|(i) 
over a range of timescales . Note that Eq. © is valid 
neither at large nor at very small z, rather at intermedi- 
ate times. 

To study the collective relaxation, we have calculated 
the collective OTCF defined as: 



cm 



{E,E/ i (e t (0)'^W)) 
<EiEi^(e*(0)-ei(0))> 



(11) 



Evidently, the calculation becomes computationally quite 
expensive even for a 1000 particle system. 

Note that the time derivative of the collective OTCF 
is directly related to the time derivative of polarizability- 
polarizability correlation function and hence, to the OKE 
signal 4]. In Fig. 0] we show the prominent power law 
regimes in the log- log plot of derivative of Cffi) at two 
temperatures just above Tin. In the case of the collective 
OTCF, the power law exponents vary only little with 
temperature. The value of the exponents are 0.33 for 
T* = 1.15 and 0.37 for T* = 1.18. 

Experimental studies on liquid crystalline systems with 
rod-like molecules find the value of the power law expo- 
nent in the range of 0.6-0.7 @. As observed before, the 
values of the exponent a are not universal. 
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FIG. 5: Power law is exhibited in the angular velocity au- 
tocorrelation function at five temperatures across I-N phase 
transition. Inset figure shows the log-log plot for three tem- 
peratures and the thick lines indicate the linear fits obtained 
in power law regimes. 
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FIG. 6: The log-log plot of the rotational mean square dis- 
placement shows the gradual onset of subdiffusive regime fol- 
lowed by a superdiffusive jump as temperature is lowered 
across I-N phase boundary. The derivative of the same plot 
has been shown in the inset figure to clearly demonstrate the 
striking behavior. At very long time the derivative should 
converge to 1 corresponding to the diffusive limit. 



B. Power law decay in angular velocity 
autocorrelation function 



We have calculated the angular velocity autocorrela- 
tion function (C^(i)) defined as: 



C u (t) 



Mo)Mt)) 

(w(0).w(0)> 



(12) 



where the angular velocity u(t) is perpendicular to the 
axis of the rotor. In Fig. [3] we have shown the behavior 
of (C u (t)) at different temperatures near Tim- We find 
that (C w (i)) has an oscillatory feature at short timescale 
due to the underdamped nature of the system. But we 
observe a pronounced power law decay with the value 
of the exponent being in the range of 1.7-1.8 at longer 
timescales. The origin of this power law can be attributed 
to the collective and correlated orientational density fluc- 
tuations. 

As in the case of orientational correlation function 
Cf (t) , a semi-quantitative understanding of (t) can be 
obtained from MCT. Since the rotational friction is given 
by Eq. © , we have the following approximate expression 
for C u {z): 



C u (z) = 



I(z + Az~ a ) 



(13) 



This expression also gives rise to the power law decay of 
C w (t) at long times. Note that decay of C u (£) occurs at 
shorter times than that of C|(t). This is the reason for 
lower amplitude of power law decay in C w (i). 



C. Rotational diffusion 

To probe the origin of the slow dynamics near the I- 
N transition, we have studied the rotational diffusion of 
the rotors on the lattice. The rotational displacement 
has been computed by integrating the ang ular velocity 
to have an unbound representation jllLI'26|. 



<t> n {t) -(f> n (0) 



dt'LU n (t') 



(14) 



In Fig. H3 the log-log plot of the mean squared angular 
displacement versus time has been shown. Interestingly, 
we observe the emergence of a subdiffussive regime at a 
timescale comparable to the plateau observed in single 
particle OTCF. This behavior becomes apparent even in 
the isotropic phase and has been illustrated by the dip 
in the derivative of the log-log plot (inset of Fig.©. The 
derivative plot shoots up suddenly after the dip and the 
motion remains superdiffussive. The system takes con- 
siderably long time to attain the diffusive limit. Our 
simulation has not been long enough to produce good 
averaging in that domain. Note that this subdiffusive 
behavior is well known in supercooled liquids and also 
has been found in calamitic systems studied with the 
Gay-Berne pair potential [Tl|. 



D. Non-Gaussian parameter 

Non-Gaussian parameters (NGP) are often useful for 
description of the dynamical heterogeneity in complex 
systems [2^, [27]]. For a system with linear rotors, the 
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FIG. 7: The rotational non-Gaussian parameter across I-N 
transition shows the increasing non-Gaussian behavior in the 
nematic phase. We don't observe any appreciable regular shift 
of the maxima with temperature as in supercooled liquids. 



rotational NGP can be denned as [Hill 0: 
<A<^) 4 } -, 



a 2 



2\2 



where 



1 N 



(15) 



(16) 



In Fig. [7] we show the time evolution of the rotational 
NGP. We observe large non-Gaussian behavior at inter- 
mediate timescales particularly at lower temperatures. 
Initially, when the motion of the rotors is ballistic, the 
rotational NGP is uniformly zero. It starts to grow at 
the comparable timescale as the appearance of the sub- 
diffusive motion and the maximum is reached when the 
rotors escape from the orientational cage and gradually 
reaches the diffusive limit at longer timescales. The for- 
mation of pseudo- nematic domains as Tjjv is approached 
from above would make the system dynamically hetero- 
geneous. The growing peak in oc^if) can, therefore, be 
ascribed to the formation of pseudo-nematic domains. 
Note that similar behavior is well known in supercooled 
liquids [Til l26l |2?1 l28l| and has been shown to be present 
in lattice models, e.g. in relaxation of phcnomenologi- 
cal Brownian rotors based on the densely frustrated XY 
model [29j . However, we do not observe any regular shift 
in the position of the maximum within the temperature 
range of our consideration. 



IV. FREE ENERGY AS A FUNCTION OF 
ORDER PARAMETER (S) 

We have computed the free energy per rotor as a func- 
tion of orientational order parameter using a variant of 



transition matrix Monte Carlo (TMMC) method recently 
proposed by Fitzgerald et al 30] . This method is consid- 
erably different than the reweighting methods for calcu- 
lation of free energy and can be incorporated in any ex- 
isting Monte Carlo simulation. An outline of the TMMC 
method applied to the present problem follows. 

The basic idea is to calculate the probability (11(5, (3)) 
that the system is in a macrostate 5 (in our case, aver- 
age orientational order parameter) at the inverse tem- 
perature (3. Note that any average observable of the 
system can be computed once the macrostate probabil- 
ity is known, e.g. in canonical ensemble, the free energy 
can be obtained as a function of order parameter from 
(3F(S,(3) = — mil (5, (3). The Boltzmann macrostate 
probability can be expressed as U(S,{3) = X^sgs 7F ( s "^)' 
where n(s,{3) is the Boltzmann probability for a partic- 
ular microstate (microscopic configuration) s. It is given 
by 7r(s,/3) = exp(— (3H S ), where H s is the value of the 
Hamiltonian for the microstate s. 

The algorithm of finding n(5, (3) is as follows: 

1. Similar to Metropolis algorithm, for a given initial 
microstate s, a new state t is proposed with prob- 
ability q Si t- As a simplification, we chose q S}t = qt, s 
though it is not strictly necessary. 

2. The probability of the move to state t being ac- 
cepted is 



r s>t (P, rj) = min{l, 



exp(rfr)Tr(t, /3) 



h 



(17) 



where 775 is the weight function corresponding to 
macrostate 5. For our purpose, we have set 775 = 
— lnn(5, f3), which corresponds to the multicanoni- 
cal approach [30j . Thus the macrostates with lower 
probability are given more weight so that toward 
the end of the simulation the histogram becomes 
flat and low probability states are sampled quite 
well. 

3. A new bookkeeping step is incorporated following 
the equilibration of the above Markov chain. At 
every step an array Cs.t (initialized to zero) is in- 
cremented as follows: 



For S^T, 
Cs.t(P) = 



Cs.t(P) 
C s ,s(P) ■ 



■r.,t(J3,V = 0) 



0)) (18) 



For S — T, C s ,s = C s ,s + 1. Note that while 
the Markov chain is guided by the multicanonical 
weight, the unweighted Boltzmann transition prob- 
abilities are stored for each visited macrostate. 

The canonical transition probability (CTP) between 
the macrostates has been calculated at some interval as 



Ps.t(P) 



Cs.t(P) 



(19) 
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FIG. 8: Free energy versus order parameter plot for a 10 x 
10 x 10 system at three different temperatures. 



The equilibrated Markov chain must obey the detailed 
balance equation Hs(P)Ps,t(P) = n T P Tj s(/3). Hence, 
the macrostate probabilities have been obtained by solv- 
ing the set of coupled linear equations iteratively. The 
weights have been updated at every 150 x L 3 steps and 
the simulations have been continued for 5 x 10 6 x L 3 steps. 
The system learned to pass through the low probability 
states automatically and reasonably good sampling was 
attained. 

The computed free energy surface is shown in Fig. [S] 
The flatness of the free energy surface near Tjn is cer- 
tainly a consequence of the very weakly first order nature 
of the transition and also partly due to the finite size 
of the system |16(. Thus, a system of this size exhibits 
large scale fluctuation in the value of the orientational 
order parameter. Such fluctuations can give rise to the 
power law decay in OTCFs 0- As the system size in- 
creases, a small barrier separating the isotropic and ne- 
matic phases appears near the critical point |lt| . This 
barrier height is directly related to the surface tension be- 
tween the isotropic and the nematic phases and it follows 
a finite size scaling law 0, |3(], E| : 



lim r-r 

L^oo 2L d - 1 



In 



n„ 



n, 



(20) 



where F s is the surface free energy, d is the dimension, L 
is the box length, and U max (P) and U min (P) correspond 
to the maximum and minimum values of macrostate 
probability respectively. 

Our result for 1000-particle system differs from one of 
the earlier reports jlfij . which uses Ferrenberg-Swendsen 
reweighting technique, as the positions of the minima 
are comparatively well-separated in our case while we 
consider temperatures far away from Tin- We feel this 
has better agreement with the average order parameter 
values that one should obtain at those temperatures. But 
for complete comparison between the methods and to 



verify the scaling relations for free energy barrier, more 
detailed study in larger systems is necessary. 

It was discussed previously that the nearly flat free en- 
ergy surface can be included in a Ginzburg-Landau type 
free energy functional |9j- One can then write down the 
following generalized Langevin equation for the fluctuat- 
ing nonconserved order parameter (SS) 0, |20| : 



d(5S) 
dt 



dtT(t-t')-^(t')+R(t), (21) 



where T is a damping coefficient, F(5S) is the Landau- 
de Gennes free energy as a function of the orientational 
order parameter and R(t) is a random velocity term re- 
lated to r by the fluctuation-dissipation theorem. As the 
temperature approaches the critical temperature T c , the 
free energy surface becomes soft. If one uses the Landau 
free energy expansion 



SF = A{T)(SS) 2 + B(T){6S) 3 + C{T)(SS) 4 



(22) 



then it can be shown that Eq.[2]can give rise to a power 
law decay of <Ag(o)> a ^ sn °rt to intermediate times. 

If the noise term in Eq. (|21|l is neglected and a Marko- 
vian approximation is made of T(t), then one finds a 
power law in the decay of (SS(0)SS(t)). This power- 
law originates from the cubic and the quartic terms in 
the free energy expansion given in Eq. I|22|l . However, 
an analytic solution in presence of the noise term R(t) 
becomes highly non-trivial and needs to be carried out 
numerically. Work in this direction is in progress. 

There are, however, two limits where one can obtain 
semi-quantitative answer directly. If the free energy sur- 
face is nearly flat, as shown in Fig. then decay of 
(6S(0)SS(i)) occurs via a generalized diffusion. In this 
case, the power law decay originates from power law be- 
havior (in time or frequency) of T(t). Note that flatness 
of free energy surface implies that A(T) m near T c . A 
study of relaxation by diffusion with a non-Markovian 
memory function was reported by Denny et al [32S |. The 
alternate limit is where a significant barrier develops be- 
tween the isotropic and the nematic phases. This, how- 
ever, is expected to give largely an exponential decay. 

Note that the Landau-de Gennes behavior is obtained 
by keeping only the quadratic term in Eq. (|22|l and the 
relaxation time is given by (in present notation): 



1 



TLdG 



2T(z = 0)A(T) 



(23) 



The exponential relaxation observed both in experiments 
and simulations follow the Landau-de Gennes behavior 
with time constant given by Eq. (|23|l . 



V. CONCLUSION 

This report contains to the best of our knowledge the 
first detailed study of dynamics in LL model near the 
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isotropic-nematic transition in the context of power law 
relaxation. Despite having a very simple intermolecular 
potential and no translational degrees of freedom, the 
LL model exhibits an array of interesting dynamical fea- 
tures near the I-N phase boundary. The orientational 
relaxation slows down at intermediate timescales possi- 
bly due to the caging effect produced by neighboring sites 
and formation of pseudo-nematic domains in the isotropic 
phase very close to Tjjv- The caging causes the rotors to 
show subdiffusive behavior. But the caging is not strong 
enough and the system remains underdamped over a long 
period of time. Essentially this leads to the vanishingly 
small free energy barrier separating the isotropic and the 
nematic phases confirming the nature of the phase tran- 
sition to be very weakly first order. 

It is interesting to note the synergy in the emergence of 
the power law decay in orientational relaxation and the 
sub-diffusive behavior of the rotational diffusion. It is 
also equally interesting to find the super-diffusive behav- 
ior that follows the sub-diffusive behavior. This signifies 
the long flights of the rotors just after they escape from 
the cage. 

The analogy between the relaxation in supercooled 
liquids and that in the isotropic phase near the I-N 
phase transition has been discussed only in recent lit- 
erature 0, El H3 ■ Even the present simple model shows 
dynamic signatures which are well-known in supercooled 
liquids. Prominent among them, are the power law re- 
laxation and the sub-diffusive behavior. It is interest- 
ing to compare with the scenario in supercooled liquid 
where one often finds two power laws, one leading to- 
wards the plateau and the second one (Von Schweidler 
law) at longer times moving away from the plateau 0. 
However, the origin of such analogous behavior in the two 
systems can be quite different. In the case of liquid crys- 
tals, the onset of long range orientational correlations and 
the formation of pseudo-nematic domains are responsible 
for the power law. Anomalies in the relaxation behavior 



in supercooled liquids, on the other hand, are believed to 
have kinetic origin. 

There now appear to exist two somewhat different in- 
terpretations of the power law decay. The first one, pro- 
posed by Gottke et al 0, [f| , is in terms of diverging ori- 
entational pair correlation function giving rise to a power 
law divergence of the memory function at low frequency 
(Eq. ©). The alternative explanation offered recently by 
Bagchi and coworkers , invokes large scale fluctuations 
in collective orientational density. The first explanation 
(in the spirit of MCT) does not require such large scale 
density fluctuation. It does, however, invoke diverging 
correlation length. The second explanation is particu- 
larly relevant for weakly first order transitions with sec- 
ond order characteristics where the two phases are sepa- 
rated by low barrier. 

Because of the length of the MD simulations required 
to obtain a reliable, statistically significant power law de- 
cay, we have been limited to study a 1000-particle system. 
The underdamped nature of the relaxation has made the 
statistical averaging demanding. However, it would be 
worthwhile to consider larger systems. We have made 
preliminary study of a (20 x 20 x 20) lattice system (8000 
particles) . We find that the relaxation behavior does not 
change significantly, but the free energy surface starts 
showing the formation of a noticeable (but still small) 
barrier at the intermediate values of the order parameter 
as expected [l6|. 
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